/*
 * =====================================================================================
 *
 *       Filename:  symmeig.h
 *
 *    Description:  Real Symmetric Matrices eigenvalue and eigenvector
 *
 *        Version:  1.0
 *        Created:  2009年06月17日 09时10分05秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Ying Wang (WY), ywang@nlpr.ia.ac.cn
 *        Company:  Institute of Automation, Chinese Academy of Sciences
 *
 * =====================================================================================
 */

#ifndef SYMMEIG_H
#define SYMMEIG_H
#include <cmath>
#include <limits>
#include "ncmatrix.h"
#include "ncvector.h"
#include "mathutils.h"

void eigsrt(NCvector<double> &d,NCmatrix<double> *v=NULL);


/**
 * \brief real symmetric matrix eigenvalue and eigvector
 *    useage: Symmeig eig(matrix); eigvalue = eig.d; eigvector=eig.z
 *
 */
struct Symmeig
{
	int n;
	NCmatrix<double> z;//the eigvector
	NCvector<double> d,e; // d the eigvalue
	bool yesvecs; // if true computer eigvector, otherwise only eigenvalue

	Symmeig(const NCmatrix<double> &a, bool yesvec=true)
		: n(a.row())
		, z(a)
		, d(n)
		, e(n)
		, yesvecs(yesvec)
	{
		tred2(); // Reduction to tridiagonal form;
		tqli();  // Eigensystem of tridiagoal matrix;
		sort();  // Sort eigvalue and eigvector by descend
	}
	Symmeig(const NCvector<double> &dd,const NCvector<double> &ee,bool yesvec=true)
	: n(dd.size())
	, z(n,n,0.0)
	, d(dd)
	, e(ee)
	, yesvecs(yesvec)
	{
	    for(int i=0;i<n;i++)
	    {
	        z[i][i] = 1.0;
	    }
	    tqli();
	    sort();
	}
	void sort()
	{
		if (yesvecs)
			eigsrt(d,&z);
		else
			eigsrt(d);
	}
	void tred2();
	void tqli();
	double pythag(const double a, const double b);

};
#endif
